#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
use Getopt::Long;
# TpA/TA
use constant DOUBLET_LENGTH => 2;

my @index_type_names = qw(TpA_TA CpA-TpG_ApC-GpT GpC_CpG CRI);
my $window = 500;
my $stepsize = $window;
my $format = 'PID';
my $index_type = 3;
GetOptions(
	   'w|window:i' => \$window,
	   's|stepsize:i'=> \$stepsize,
	   'f|format:s' => \$format,
	   'i|indextype:i' => \$index_type,
	  );

if( $format =~ /^WIG/i ) {
    $format = 'WIG';
} elsif( $format =~ /PID/i || $format =~ /percent/ ) {
    $format = 'PID';
} else {
    die("Don't understand format $format\n");
}

my @doublets;# = qw(CA TG AC GT CG GC);

for my $base1 ( qw(C A G T) ) {
    for my $base2 ( qw( C A G T) ) {
	push @doublets, $base1.$base2;
    }
}

for my $file ( @ARGV ) {
    my $in = Bio::SeqIO->new(-format => 'fasta',
			     -file   => $file);
    # drop some suffixes
    $file =~ s/\.(fa|fasta|seq|dna|nt)$//;
    my $ofh;
    if( $format eq 'WIG' ) {
	open($ofh => ">$file.RIP_".$index_type_names[$index_type].".wig") || die $!;
	printf $ofh "type track=%s name=\"%s\" description=\"%s %s\"\n",
	'wiggle_0', 'RIP index',
	$index_type_names[$index_type], 
	'RIP index calculation';
    } else {	
	open($ofh => ">$file.RIP.dat") || die $!;
	print $ofh join("\t", qw(CHROM WINDOW), @doublets, @index_type_names),"\n";
    } 
    while( my $seq = $in->next_seq ) {
      if ( $seq->length < $window ) { 
	# skip sequences that are shorter than the window size
	next;
      }
	my $chrom = $seq->display_id;
	my $seqstr = $seq->seq;
	my $i = 0;
	my $j = 0;
	my $len = $seq->length;
	if( $format eq 'WIG' ) {
	    printf $ofh "fixedStep chrom=%s start=%d step=%d\n",
	    $chrom, 1, $stepsize;
	}
	while( $i <= $len ) {
	    my %freq;
	    for my $d ( @doublets ) { 
		$freq{$d} = 0;
	    }
	    my $r;
	    if( $i > ( $len - $window) ){ 
	      # skip the end of the contig
	      #next;
	      # do some wrap around for the end-case
	      # -------------------->
	      #          |---------|
	      # |---------|
	      $r = substr($seqstr, $i-$window, $window);
	    } else {
	      $r = substr($seqstr,$i, $window);
	    }
#	    my $TAn = ( $r =~ tr/[ATat]/[ATat]/);
#	    my ($TA) = 0;
	    my $wlen = length($r);
	    for( my $k = 0; $k < $wlen - DOUBLET_LENGTH; $k++) {
		my $w = uc(substr($r,$k,DOUBLET_LENGTH));
		$freq{$w}++;
#		$TA +=2 if ($w eq 'TA' || $w eq 'AT');
	    }
	    $i+= $stepsize;
	    my @rip_indexes;
	    push @rip_indexes, ($freq{AT}) ? $freq{TA}/$freq{AT} : -1;
	    push @rip_indexes, ($freq{AC} + $freq{GT}) ?
		($freq{CA} + $freq{TG})/($freq{AC} + $freq{GT}) : -1;
	    push @rip_indexes, $freq{CG} ? ($freq{GC} / $freq{CG}) : -1;
	    # CRI is substrate - product;
	    push @rip_indexes, $rip_indexes[1] - $rip_indexes[0];
	    if( $format eq 'PID' ) {
		print $ofh join("\t", $chrom,$j++, 
				(map { $freq{$_} } @doublets),
				map { sprintf("%.2f",$_) } @rip_indexes),"\n";
	    } elsif( $format =~ /WIG/ ) {
		printf $ofh "%.2f\n",$rip_indexes[$index_type];
	    }
	}
    }
}

